Stability of Discrete Solitons in the Presence of Parametric Driving 
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In this brief report, we consider parametrically driven bright solitons in the vicinity of the anti- 
continuum limit. We illustrate the mechanism through which these solitons become unstable due 
to the collision of the phase mode with the continuous spectrum, or eigenvelues bifurcating thereof. 
We show how this mechanism typically leads to complete destruction of the bright solitary wave. 
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Introduction. In the past few years, differential- 
difference dispersive equations where the evolution vari- 
able is continuum but the spatial variables are discrete, 
have been the focus of intense research efforts jf,]. The 
key reason for the increasing interest in this research di- 
rection can be attributed to the wide range of pertinent 
applications ranging, from e.g., the spatial dynamics of 
optical beams in coupled waveguide arrays in nonlinear 
optics , to the temporal evolution of Bose- Einstein con- 
densates (BECs) in deep, optically-induced, lattice po- 
tentials in soft-condensed matter physics 0, or even to 
the DNA double strand in biophysics |4| among others. 

One of the key models that has emerged in all of the 
above settings, either as describing e.g., the envelope 
wave of the electric field in the optical setting p|, or 
describing the wavefunction at the nodes of the optical 
lattice in BECs 6] , is the discrete nonlinear Schrodinger 
(DNLS) equation. This prototypical lattice model fea- 
tures a dispersive coupling between nearest-neighbors, 
and a cubic onsite nonlinearity. 

The above spatially discrete model bears a number 
of interesting similarities and differences, in comparison 
with its continuum sibling, the famous (integrable in 1- 
spatial dimension) nonlinear Schrodinger equation (NLS) 
• One of the key differences is the breaking of one of the 
important invariances of the NLS model, namely of the 
translational invariance that is responsible for momen- 
tum conservation in that setting. On the contrary, the 
discrete model carries an integer-shift invariance. This 
has some important implications, among other things, to 
the nature of the solutions of the discrete model. In fact, 
it was realized through perturbative calculations y| and 
subsequently more rigorously justified 9] that the princi- 
pal (single-humped solitary wave) solutions of the latter 
model can only be centered on a lattice site or between 
two lattice sites. In the continuum case, the center of 
the solution is a free parameter due to the continuum 



invariance. 



On the other hand, one of the important similarities of 
the discrete model to the continuum one is the presence 
of the so-called phase or gauge invariance (which is asso- 
ciated with the overall freedom of selecting the solution's 
phase). The conservation law related to this invariance 
is the one of the L 2 (respectively I 2 ) norm, or "mass" of 
the solution. This invariance is the main focal point of 
the present work. In particular, we introduce, arguably, 
the simplest possible perturbation that breaks the rele- 



vant invariance, in the form of a parametric drive. The 
relevance of such a term involving a perturbation pro- 
portional to the complex conjugate of the field has been 
discussed in a variety of earlier works (see e.g. |lf)j | and 
references therein) . A specific physical setting where this 
type of perturbation arises can be found by looking at the 
envelope equation of a system of parametrically driven 
(undamped) coupled torsion pendula as discussed in 
(with the difference that the envelope wave expansion 
should be performed in a genuinely discrete setting sim- 
ilarly to [12j rather than near the continuum limit as in 
[Tlj)- The aim of this exposition is to examine how the 
breaking of this invariance results in an eigenvalue that 
bifurcates from the origin of the spectral plane, when 
linearizing around the most fundamental, solitary wave 
solution. We argue (analytically and support numeri- 
cally) that this eigenvalue can lead to an instability of 
the solitary wave for an isolated value of the parametric 
drive even at the so-called anti-continuum limit where 
lattice sites are uncoupled. For non-vanishing couplings, 
the same eigenvalue leads to a wide interval of paramet- 
ric instabilities in the two-parameter space (of parametric 
drive versus inter-site coupling) that we explore both an- 
alytically and numerically. Within this interval, we also 
elucidate the typical numerical behavior of the solitary 
wave solutions, using direct numerical simulations of rel- 
evant unstable waveforms. 

Our presentation will be structured as follows. In the 
next section, we present our analytical setup and pertur- 
bative results. Then, we compare our analytical findings 
with the results of numerical computations. Finally, we 
summarize our findings and present our conclusions, as 
well as motivate some questions for future study. 

Setup and perturbation analysis. The model we con- 
sider is the perturbed (i.e., parametrically driven) dis- 
crete non-linear Schrodinger equation of the form 



i(j> n = -CA 2 <j) n 



J 2 </>n + A0„ +70„, 



(1) 



where C is the coupling constant between two adjacent 
sites of the lattice, A 2 0„ = {(f> n+ i — 2(f> n +(f> n ^i) is the dis- 
crete Laplacian, A is the propagation constant in optics 
or the chemical potential in BECs, and 7 is the strength 
of the parametric drive. 

We focus our attention on a standing wave profile so 
that <p n is time- independent. In this case, <p n satisfies 



-CA 2 



<pl + A4> n + 70„ = 0. 



(2) 



2 



In the uncoupled (or so-called anti-continuum) limit of 
C = 0, the solution of is (f) n = 0, ±yA + 7. We exam- 
ine here the most fundamental single-hump solitary wave 
solutions which in the anti-continuum limit emanate from 
a single-site excitation of the form: 



= 0, n ^ 0, u° = VA + 7- 



(3) * 



The continuation of J2J for small coupling C can be 
calculated analytically through a perturbative expansion. 
By substituting into the steady state equation (J2J) u n = 
u n + Cu^ + C 2 !^ + . . . , one can calculate that up to order 
0(C 2 ) 



v / a~ T7 + c/yi+7, n = o, 

»,, = < C7VA + 7, n = -1,1, (4) 

0, 72 otherwise. 

To perform linear stability analysis to the discrete soli- 
tary waves of the form of Eq. |0J , we introduce the fol- 
lowing linearization ansatz 



Se r , 



Substituting into Q yields the following linearized equa- 
tion to 0{8) 



= -CA 2 e n - 21 



u 2 n e n + Ae„ + 7e„ 



(5) 



Writing e n (t) = r\ n +i£n and assuming that u n is real, 
eq. gives (see, e.g., |l3||') 



'In 

L 



£+(C) 
-£-(<?) 



H 



,(6) 



where the operator £_(C) and £+(C) arc defined as 
£_(C) = -CA 2 - (3u* - A - 7) and £+(C) = -C*A 2 - 
(m„ — A + 7) . The stability of w n is then determined by 
the eigenvalues of H. 

Let the eigenvalues of Ti. be denoted by ioj, which im- 
plies that u n is stable if Im(w) = 0. Because (0 is linear, 
we can eliminate one of the 'eigenvectors', for instance £ n , 
from which we obtain the following eigenvalue problem 



£ + (C)£_(C)t7„ = u> 2 r) n = Vt-q n . 



(7) 



As before, we expand the eigenvector rj n and the eigen- 
value f2 as 

r) n = r?° + Cri + 0(C 2 ), n = n° + cn l + 0(C 2 ). 

Substituting into Eq. (JJJ) and identifying coefficients 
for consecutive powers of the small parameter C yields 



[£ + (0)£_(0)-fi°K = 0, 
[£ + (0)£_(0)-fi ]^ = /, 



with 

/=[-(A 2 



2«)£_(0)-£ + (0)(A 2 + 6u°« 1, 



(8) 
(9) 



■ ^] r,° n 
(10) 




FIG. 1: The smallest eigenvalue for two values of 7, namely 
7 = 0.02 and 7 = 0.1. The dashed lines are the approximate 
analytical estimate of the relevant frequency from Eq. 11121 . 
The lower curves correspond to 7 = 0.02. 



First, let us consider the order 0(1) equation (J8J. One 
can do a simple analysis to show that there are only two 
eigenvalues, i.e., Q° = A 2 — 7 2 and f2° = 4(A + 7)7. 
f2° = A 2 — 7 2 has infinite multiplicity and is related 
to the continuous spectrum that will be discussed later. 
Therefore, our interest is in Sl° = 4(A + 7)7 that has 
the normalized eigenvector 77^ = 0, n 7^ and t/q = 1. 
This eigenvalue is the formerly zero eigenvalue due to 
the phase or gauge invariance of the DNLS equation in 
the absence of parametric driving. 

The continuation of the eigenvalue il° = 4(A + 7)7 
when the coupling C is turned on can be calculated from 
Eq. ©. Due to the corresponding eigenvector having 
?7„ = for n 7^ 0, we only need to consider the site 
n = 0. In this case, / = —87 + SI 1 . The solvability 
condition of Eq. © using, e.g., the Fredholm alternative 
requires / — from which we immeadiately obtain that 
SI 1 = 87. Hence, the smallest eigenvalue of a one-site 
discrete soliton solution of Eq. Q is 



Q = 4(A + 7 )7 + 8 7 C7 + 0(C7 2 ), 



(11) 



or 



= ±2v/(A + 7 ) 7 ±2- 



If 7 = 0(C), Eq. becomes 



C + 0(C 2 ). (12) 



u = ±2 V /A^v / C + 0(C). 



(13) 



Next, we have to proceed with calculating the contin- 
uous spectrum of the operator £ + (C)£_(C) J7J. When 
C = 0, all the continuous spectrum of the operator lies 
at SI = A 2 — 7 2 as was mentioned before. When C is 
increased, the eigenvalues spread on the imaginary axis 
creating a phonon band. Using a plane wave expansion 
rj n = ae lKn + be~ tKn yields the dispersion relation 

SI = (A + 7 + 2C - 2Ccosk)(A - 7 + 2C - 2Ccosk). 
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FIG. 2: The stability-instability region in the two-parameter 
space 7 — C. The solid lines give the numerically obtained 
separatrices, while the dash-dotted and dashed ones the ana- 
lytical approximations of Eqs. 114H and 1151 respectively. 
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FIG. 3: The eigenvalue structure of a single-hump solitary 
wave for 7 = 0.02 and C = 1.0 (top left panel), as well as C = 
1.4 (top right panel). The bottom panel shows the trajectory 
of one of the unstable eigenvalues as C changes. 

Hence, the continuous band lies between £1^ = A 2 — j 2 
(when k = 0) and flu = A 2 - 7 2 + 8C(A + 2C) (when 

K = 71"). 

For small 7, the instability of a one-site discrete 
breather of is caused by the collision of the small- 
est eigenvalue with an eigenvalue bifurcating from 
flL- However, here we assume that the bifurcating eigen- 
value does not move very fast in the spectral plane such 
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FIG. 4: (Color online) The spatio-temporal evolution of an 
unstable single-hump solitary wave at 7 = 0.02 and C = 1.4. 
The contour plot of square modulus |0„| 2 is shown. 



that it can be represented by flL- For large 7, the in- 
stability is due to the collision of the smallest eigenvalue 
and f^L- Equating those quantities will give the critical 
7 as a function of the coupling constant C, i.e. 

7c 1 , - -^A - + i V9A 2 + 16C(A + C), (14) 

7 2 r = --A --C+ - ^9A 2 + 56CA + 96C 2 . (15) 
5 5 5 

The two approximate 7* r above coincide at C = and 
7* = (V2- l)/2 w 0.2071. Notice that at that level 
the relevant calculation is analytically exact (i.e., there 
is no approximation and the solitary excitation will be 
unstable for C — only for 7 = 7*. 

Numerical results. We now proceed to testing our ana- 
lytical results for the parametrically driven discrete non- 
linear Schrodinger system numerically. We start by ex- 
amining the validity of our analytical prediction for the 
eigenfrequency corresponding to the phase mode which 
bifurcates from li — because of the presence of the para- 
metric drive according to the expression (|12|) . Figure ^ 
shows this prediction as a function of C for two differ- 
ent values of 7. Clearly the prediction is fairly accurate 
for small C and its range of validity is wider for smaller 
values of 7. 

We now turn to the examination of the two parameter- 
plane of the parametric drive 7 versus the coupling 
strength C. Figure [21 provides a full description of the 
dynamics of the parametrically driven DNLS model re- 
garding the intervals of stability/instability of the most 
fundamental, single-hump solitary wave solution of the 
model. The solid lines show the numerically obtained 
separatrices between the stable and unstable parametric 
regimes of the model, while the dashed and dash-dotted 
lines give the analytical prediction for the stability range 
as obtained by the conditions of collision of the phase 
mode eigenfrequency with the continuous spectrum from 
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Eqs. ljl4"|) - (fH)|l . We observe that the prediction of Eq. 
(|15fl is in very good agreement with the numerical obser- 
vations for the occurrence of the instability point. This 
is because the collision typically occurs indeed with the 
upper band edge of the continuous spectrum (rather than 
with an eigenvalue bifurcating from it) and also typically 
the collision occurs for small C for which the analytical 
approximation of Eq. i|12[l is a very good approximation. 
On the other hand, the slightly less satisfactory agree- 
ment with the prediction of Eq. ifTil) occurs due to the 
collision with eigenvalues bifurcating from the lower edge 
of the continuous spectrum (see also Fig. [3] below) and 
also for relatively large C's for which higher order terms 
in the expansion of l(T2*|l should be expected to contribute. 

Figure illustrates the typical instability scenario for 
weak parametric drives (7 = 0.02 in this figure). As 
C increases, the eigenvalue which is associated with the 
phase mode moves towards the continuous spectrum (top 
left panel). Eventually for C ~ 1.021 it collides with an 
eigenvalue pair that has bifurcated from the lower band 
edge of the continuous spectrum. Due to the opposite 
Krein signature of these eigenvalues (see e.g. the relevant 
discussion in 0]), their collision leads to an oscillatory 
instability and the bifurcation of a complex quartet of 
eigenvalues (top right panel) . Eventually, as is shown in 
the bottom panel, the eigenfrequencies return to the real 
axis to re-stabilize the configuration for C > 1.78. 

One can also notice from Fig. [21 that there is a mini- 
mum 7 Jn below which the soliton is stable all the way to 
the continuum limit. Numerically, j m 1=3 0.0135. When 
7 is less than 7 m the eigenvalue that moves towards 
the continuous spectrum does not collide with the eigen- 
value bifurcating from the phonon band Vtr,- Instead, it 
decreases before the collision occurs. This shows that 
the second order correction 0(C 2 ) of pif! dominates the 
leading order expression. 



We now turn to the examination of the dynamical be- 
havior of the unstable solutions obtained above. The 
direct numerical evolution of an unstable solution of Eq. 
Q is shown in Figure 01 We have confirmed that this 
dynamics is typical of the unstable parameter range. The 
figure shows that eventually the solution becomes subject 
to the oscillatory instability that was illustrated in Fig. 
[3] and is ultimately destroyed completely. This may also 
be expected on the basis of the fact that this is the fun- 
damental coherent structure solution and for the same 
parameter set there appears to be no other stable dy- 
namical state (other than <fi n = 0) to which the initial 
condition may transform. 

Conclusions. In this short communication, we visited 
the topic of parametrically driven lattices of the nonlin- 
ear Schrodinger type. We have shown that the dynamics 
of these lattices is considerably different than those of 
the regular DNLS equation. This is due to the driving- 
induced bifurcation of the phase mode (associated with 
the gauge invariance of the NLS equation). Collision of 
this mode with eigenfrequencies stemming from the con- 
tinuous spectrum leads to a wide parametric regime of in- 
stabilities of the fundamental solitary wave in this model. 
Our perturbative analysis captures quite accurately the 
relevant eigenvalue (especially for weak couplings) and 
provides a fair estimate of the instability threshold in the 
parameter-space of the system. The result of the ensuing 
oscillatory instability is the destruction of the fundamen- 
tal soliton, a feature absent from the regular DNLS model 
(where this solution is stable for all parameter values). 

It would be interesting to expand the present consid- 
erations to other variants of the discrete parametrically 
driven model such as its higher-dimensional analogs, the 
defocusing case, or also damped variants of these lattice 
models. Such considerations are currently under study 
and will be reported in future publications. 
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